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Abstract 



We present a parametric deformable model which recovers image components with 
a complexity independent from the resolution of input images. The proposed model 
also automatically changes its topology and remains fully compatible with the gen- 
eral framework of deformable models. More precisely, the image space is equipped 
with a metric that expands salient image details according to their strength and 
their curvature. During the whole evolution of the model, the sampling of the con- 
tour is kept regular with respect to this metric. By this way, the vertex density is 
reduced along most parts of the curve while a high quality of shape representation 
is preserved. The complexity of the deformable model is thus improved and is no 
longer influenced by feature-preserving changes in the resolution of input images. 
Building the metric requires a prior estimation of contour curvature. It is obtained 
using a robust estimator which investigates the local variations in the orientation of 
image gradient. Experimental results on both computer generated and biomedical 
images are presented to illustrate the advantages of our approach. 

Key words: deformable model, topology adaptation, resolution adaptation, 
curvature estimation, segmentation/reconstruction. 
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1 Introduction 



In the field of image analysis, recovering image components is a difficult task. 
This turns out to be even more challenging when objects exhibit large varia- 
tions of their shape and topology. Deformable models that are able to handle 
that kind of situations can use only little a priori knowledge concerning im- 
age components. This generally implies prohibitive computational costs (see 
Table OP). 

In the framework of parametric deformable models, most authors pll2|l3] pro- 
pose to investigate the intersections of the deformable model with a grid that 
covers the image space. Special configurations of these intersections character- 
ize the self collisions of the mesh. Once self-instersections have been detected, 
local reconfigurations are performed to adapt the topology of the model ac- 
cording to its geometry. To take advantage of all image details, the grid and 
the image should have the same resolution. An other method consists in 
constraining the lengths of the edges of the model between two bounds. Self- 
collisions are then detected when distances between non-neighbor vertices fall 
under a given threshold. Topological consistency is recovered using local op- 
erators that reconnect vertices consistently. Using all image details requires 
edges to have the same size as image pixels. The complexities of all these 
methods are thus directly determined by the size of input data. 

In the framework of level-set methods, boundaries of objects are implicitly 
represented as the zero level set of a function / [5ll6H71l8 ]. Usually / is sampled 
over a regular grid that has the same resolution as the input image. Then 
/ is iteratively updated to make its zero level-set approach image contours. 
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Table 1 

Complexities of the deformable models that automatically adapt their topology. The 
length of the deformable model is denoted I, the with of a pixel is denoted d. The 
size of the image is denoted |/|, and k denotes the width of the narrow band when 
this optimization is used. This shows that the complexities of these algorithms are 
completely determined by the resolution of the input image. 



Model 


Complexity per iteration 


T-Snake[3] 




Simplex mesh[l] 


0(h) 


Distance constraints [^) 


Oaiog^)) 


Level-set |5|7] 




Level-set with narrow band [9j 


O(fc^)+0(^logQ)) 



Even with optimization methods which reduce computations to a narrow band 
around evolving boundaries PITU] . the complexity of these methods is deter- 
mined by the resolution of the grid and hence by the resolution of the input 
image. 

In [TT] a method is proposed to adapt the resolution of a deformable model 
depending on its position and orientation in the image. The main idea is to 
equip the image space with a Riemannian metric that geometrically expands 
parts of the image with interesting features. During the whole evolution of 
the model, the length of edges is kept as uniform as possible with this new 
metric. As a consequence, a well chosen metric results in an accuracy of the 
segmentation process more adapted to the processed image. 

In this first attempt the metric had to be manually given by a user. This 
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was time consuming and the user had to learn how to provide an appropriate 
metric. Our contribution is to propose an automated way of building metrics 
directly from images. The accuracy of the reconstruction is determined by the 
geometry of image components. The metric is built from the image 

(1) to optimize the number of vertices on the final mesh, thus enhancing the 
shape representation, 

(2) and to reduce both the number and the cost of the iterations required to 
reach the rest position. 

Property ([1]) is obtained by building the metric in such a way that the length 
of the edges of the model linearly increases with both the strength and the 
radius of curvature of the underlying contours. 

Property ([2]) is ensured by designing the metric in such a way that a coarse 
sampling of the curve is kept far away from image details, while it progressively 
refines when approaching these features. 

To build a metric which satisfies these constraints, the user is asked for only 
three parameters: 

• Sref: the norm of the image gradient over which a contour is considered as 
reliable, 

• ^max^ the maximum length of the edges of the deformable model (this is 
required to prevent edges from growing too much which would lead to nu- 
merical instability), 

• /jnin: the minimum length of an edge (typically this corresponds to the half 
width of a pixel) . 

Over a sufficient resolution of the input image, the gradient magnitude as well 
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as the curvature of objects do not depend on the samphng rate of the input 
image. Consequently the computational complexity of the segmentation algo- 
rithm is determined only by the geometrical complexity of image components 
and no longer by the size of input data. 

The dynamics of the deformable model is enhanced too. In places without 
image structures the length of edges reaches its maximum. This results in 
(i) less vertices and (ii) larger displacements of these vertices. By this way 
both the cost per iteration and the number of iterations required to reach 
convergence are reduced. 

We point out that these enhancements do not prevent the model from changing 
its topology and that the complexity of these topology changes is determined 
by the (reduced) number of vertices. Other methods, such as those presented 
in [1] and [3] require the model to be resampled (at least temporarily) on 
the image grid. As a result, these methods cannot take advantage of better 
sampling of the deformable curve to reduce their complexity. 

Fig. [T] illustrates the main idea of our approach and offers a comparison with 
the classical snake approach and with a coarse-to-fine snake method. Note 
that the same parameters are used for all three experiments: force coefficients, 
initialization, convergence criterion. First, it appears clearly that our approach 
achieves the same segmentation quality as regular snakes with a high preci- 
sion. A coarse-to-fine approach may fail to recover small components. Second, 
computation times are greatly improved with our approach (about 6 times 
faster). The coarse-to-fine approach is also rather slow since a lot of time is 
spent extracting details that were not present at a coarser level. Third, the 
number of vertices is optimized according to the geometry of the extracted 
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shape (3 times less vertices). 



Moreover, the proposed model remains compatible with the general frame- 
work of deformable models: to enhance the behavior of the active contour, 
any additional force p^l3lll4] may be used without change. In practice, with 
the same sets of forces, the visual quality of the segmentation is better with an 
adaptive vertex density than in the uniform case. Indeed, along straight parts 
of image components, the slight irregularities that result from the noise in the 
input images are naturally rubbed out when fewer vertices are used to repre- 
sent the shape. In the vicinity of fine image details an equivalent segmentation 
accuracy is achieved in both the adaptive and uniform cases. 



At last, the approach proposed to build the metric is almost fully automated. 
Thus, only little user interaction is required. However it remains easy to in- 
corporate additional expert knowledges to specify which parts of image com- 
ponents have to be recovered accurately. 



This paper is structured as follows: in Sect. [2] we describe the proposed de- 
formable model and we show how changing metrics induces adaptive resolu- 
tion. In Sect. [31 we explain how suitable metrics are built directly from images 
using a robust curvature estimator. Experimental evaluation of both the pro- 
posed model and the curvature estimator are presented in Sect. HI 



7 



2 Deformable Model 



2.1 General description 

Our proposed deformable model follows the classical energy formulation of 
active contours [15j: it is the discretization of a curve that is emdedded in the 
image space. Each of its vertices undergoes forces that regularize the shape 
of the curve, attract them toward image features and possibly tailor the be- 
havior of the model [T2p3] for more specific purposes. In this paper, classical 
parametric snakes are extended with the ability to (i) dynamically and au- 
tomatically change their topology in accordance with their geometry and (ii) 
adapt their resolution to take account of the geometrical complexity of re- 
covered image components. In addition, only little work is necessary to adapt 
the proposed model for three-dimensional image segmentation (see [16] for 
details). In this three dimensional context, the number of vertices is further 
reduced and computation times are consequently greatly improved. 

2.2 Resolution adaptation 

During the evolution of the model, the vertex density along the curve is kept 
as regular as possible by constraining the length of the edges of the model 
between two bounds 5 and C,5: 

S<Le{u,v)<(S , (1) 

In Le denotes the length of the line segment that joins u and v. The 
parameter 6 determines lengths of edges and hence vertex density along the 
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curve. The parameter ( determines the allowed ratio between maximum and 
minimum edge lengths. 



At every step of the evolution of the model, each edge is checked. If its length 
is found to be less than 6 then it is contracted. In contrast, if its length exceeds 
the (6 threshold then the investigated edge gets split. To ensure the conver- 
gence of this algorithm ( must be chosen greater than two. In the following 
( is set to the value 2.5, which provides satisfying results in practice. The 
parameter S is derived from the maximum edge length /max specified by the 
user as 5 = /max/C- 

Adaptive resolution is achieved by replacing the Euclidean length estimator 
by a position and orientation dependent length estimator L/j in ([1]). In 
places where underestimates distances, estimated lengths of edges tend to 
fall under the 6 threshold. As a consequence, edges tend to contract and the 
resolution of the model locally decreases. In contrast, the resolution of the 
model increases in regions where overestimates distances. 

More formally, Riemannian geometry provides us with theoretical tools to 
build such a distance estimator. In this framework, the length of an elementary 
displacement ds that starts from point {x, y) is expressed as: 



where G associates a positive-definite symmetrical matrix with each point of 
the space. The G mapping is called a Riemannian metric. From ([2]) follow the 
definitions of the Riemannian length of a path as 



ds||^ = *ds X G{x,y) x ds. 



(2) 




(3) 
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and of the Riemannian distance between two points u and v as 

dR{u,v) = inf Lij(7) , (4) 

76C 

where C contains all the paths that join u and v. It is thus easily seen that 
defining the G mapping is enough to completely define our new length esti- 
mator L/j. How this mapping is built from images to enhance and speed up 
shape recovery is discussed in Sect. |3l 

2.3 Topology adaptation 

During the evolution of the model, care must be taken to ensure that its 
interior and exterior are always well defined: self-collisions are detected and 
the topology of the model is updated accordingly (see [TT] for more details on 
topology adaptation). 

Since all the edges have their length lower than (6, a vertex that crosses over 
an edge {u,v) must approach either u or v closer than ^{(6 + dmax), where 
dmax is the largest distance covered by a vertex during one iteration. Self- 
intersections are thus detected by looking for pairs of non-neighbor vertices 
{u, v) for which 

dE{u, V) < ^(C*^ + dmax) ■ (5) 

It is easily shown that this detection algorithm remains valid when c?^; is 
replaced with a dji distance estimator as described in Sect. 12. 2[ With a naive 
implementation, the complexity of this method is quadratic. However, it is 
reduced to 0(?2logn) by storing vertex positions in an appropriate quadtree 
structure. 

Detected self-intersections are solved using local operators that restore a con- 
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sistent topology of the mesh by properly reconnecting the parts of the curve 
involved in the collision. 



2.4 Dynamics 

Theoretically, in a space equipped with a Riemannian metric, the position x 
of a vertex that undergoes a force F follows equation 

mxk = Fk-Y^ T^jXiXj , (6) 

where the Ff^ coefficients are known as the Christoffel's symbols associated 
with the metric: 

2 V V-^^;, dx, dxi) ■ ^ ' 

However, the last term of ([6]) is quadratic in x and has therefore only little 
influence when the model is evolving. Furthermore, once at rest position it 
cancels and has consequently no impact on the final shape. Therefore it is 
neglected and we get back the classical Newton's laws of motions. Experi- 
mentally, removing this term does not induce any noticeable change in the 
behavior of the model. 



3 Tailoring Metrics to Images 

3.1 Geometrical interpretation 

For any location (x, y) in the image space, the metric G{x, y) is a positive- 
definite symmetrical matrix. Thus, in an orthonormal (for the Euclidean norm) 
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base (vi,V2) of eigenvectors, G{x,y) is diagonal with coefficients (/ii,/i2). 
Hence, the length of an elementary displacement ds = xiVi + X2V2 is ex- 
pressed as 

||ds||^ = nixl + fi2xl . (8) 
This shows that changing the Euclidean metric with a Riemannian metric 
locally expands or contracts the space along vi and V2 with ratios 1 / ^/JIi and 
Suppose now that Le replaced by Lji in ([T]). In a place where the 
edges of the model are aligned with vi this yields 

^ <LE{e)<^. (9) 



//ii ^/ii 

Of course a similar inequality holds in the orthogonal direction. This shows 
that (from a Euclidean point of view) the vertex density on the mesh of the 
model is increased by a ratio in the direction of vi and by a ratio y/Jl2 
in the direction of V2. Therefore, at a given point of the image space, a direct 
control over the vertex density on the deformable mesh is obtained by properly 
tweaking vi, V2, /ii and ^2 in accordance with underlying image features. 

Although these eigenvectors and eigenvalues could be given by a user, it is a 
tedious and complicated task. It is therefore much more attractive and efficient 
to have them selected automatically. The subsequent paragraphs discuss this 
problem and describe a method to build the metric directly from the input 
image in such a way that the vertex density of the mesh adapts to the geometry 
of image components and no longer depends on the resolution of input data. 

This property is interesting because the model complexity is made indepen- 
dent from the image resolution and is defined instead only by the geometrical 
complexity of the object to recover. Now, the geometrical complexity of an ob- 
ject embedded in an image cannot exceed the image resolution. Furthermore, 



12 



since objects do not have high curvatures everywhere on their boundary, this 
complexity is generally much smaller. 

3.2 Definition of metrics 

Two cases have to be considered: 

(1) the case for which the model has converged, and for which we expect it 
to follow image contours, 

(2) and the case for which the model is still evolving. Thus parts of the curve 
may be far away of image details or cross over significant contours. 

In casein the length of its edges is determined by (i) the geometrical properties 
of the recovered image components and (ii) the certainty level of the model 
position. More precisely, the length of edges is an increasing function of both 
the strength and curvature of the underlying contours. 

In case |2l two additional sub-cases are possible. 

• If the model crosses over the boundary of image components, the vertex 
density on the curve is increased. By this way the model is given more 
degrees of freedom to get aligned with the contour. 

• In a place with no image feature [i.e. far away from image contours) vertex 
density is kept as low as possible. As a consequence, the number of vertices 
and hence the computational complexity decreases. Moreover since edge 
length is increased, vertices are allowed to travel faster. The number of 
iterations required to reach the rest position of the model is thus reduced. 
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To obtain these properties, the eigenstructure of the metric is chosen as follows 



Vi = n and fii 
V2 = n-*- and fi2 



2 ^ J" 

■Srcf ^ref 



-I 1,- 



(10) 



X fii 



where n denotes a vector normal to the image contour, and the notation [■]„ ;, 
constrains its arguments between the bounds a and b. The parameters s and k, 
respectively denote the strength and the curvature of contours at the investi- 
gated point of the image. The parameter Kmax corresponds to the maximum 
curvature that is detected in the input image. The different possible situations 
are depicted on Fig. [2l Computing these parameters directly from images is 
not straightforward and is discussed in Sect. 13.41 

The parameter Sref is user-given. It specifies the strength over which a contour 
is assumed to be reliable. If an edge runs along a reliable contour, then its 
length is determined only by the curvature of the contour (see region B in 
Fig. [2] and Fig. [3]- left). In other cases the length of edges decreases as the 
contour gets weaker. 

The parameter reference curvature for which the length of edges 

is allowed to vary between 6 and (6 only. Below this curvature contours are 
assumed to be straight and the length of edges remains bounded between 6 and 
(6. Along more curved contours, the length of edges increases linearly with 
the radius of curvature (see Fig. [3]-left). This parameter is easily computed 
from the minimal length allowed by the user for the edges (see Fig. [2]): 



^ref 



S = L 
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where Iram denotes the chosen minimal length. To take advantage of all the 
details available in the image, it is usual to set /min to the half width of a pixel. 

Note that all the parameters are squared to compensate for the square root 
in (Q. By this way the length of edges varies linearly with both 1/s and 1/k. 

3.3 Influence of the resolution of input images 

For input images with sufficiently high sampling rates, both s and k are deter- 
mined only by the geometry of image components. Consequently, the vertex 
density during the evolution of the model and after convergence are completely 
independent from the image resolution (see experimental results on Fig. [TO!) . 

If the sampling rate is too low to preserve all the frequencies of objects, con- 
tours are smoothed and fine details may be damaged. As a result, s is under- 
estimated over the whole image and k along highly curved parts of objects. In 
these areas, the insufficient sampling rate induces longer edges and details of 
objects cannot be represented accurately. However, these fine structures are 
not represented in input images. As a consequence, it is not worth increasing 
the vertex density since small features cannot be recovered even with shorter 
edges. 

In featureless regions, or along straight object boundaries, the length of the 
edges depend neither on s nor on k and remains bounded between 6 and (6 
(see Fig. 12]). As a result the improper sampling rate of the image has no impact 
on the vertex density on the deformable curve in these regions, which remains 
coarse. 
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As a result of this behavior, the segmentation process is able to take advantage 
of all finest details that can be recovered in the image. Indeed, when the 
resolution of the input image is progressively increased, contours and details 
are restored and s and k get back their actual value. As a consequence the 
lengths of edges progressively decrease in image parts with fine details while 
remaining unchanged elsewhere. At the same time the global complexity of 
the model increases only slightly with the resolution of images until all the 
frequencies of image components are recovered. If the image is oversampled, s 
and K are left invariant, and the number of vertices, and hence the complexity 
remains unaffected. 

These properties are illustrated experimentally on Fig. [TOl Fig. [11] and [T2j 

3.4 Computing strength and curvature of contours from images 

To tailor metrics to enhance and fasten image segmentation we need to esti- 
mate both the strength s of image contours and their curvature k. 

Consider a unit vector v and Qv{x, y) = {v ■ VI(a;, y j)^- This quantity reaches 
its maximum when v has the same orientation (modulo vr) as VI{x,y). The 
minimum is reached in the orthogonal direction. To study the local variations 
of the image gradient it is convenient to consider the average of over a 
neighborhood. It is expressed in a matrix form as 

g^(x, = V X VI X *VI X V , (12) 

where (■) denotes the average of its argument over a neigborhood of point 
{x,y). The positive-definite symmetrical matrix J = VI x *VI is known as 
the gradient structure tensor. This operator is classically used to analyze lo- 
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cal structures of images [T7], since it characterizes their local orientations. It 
is further used for texture and image enhancement in anisotropic diffusion 
schemes [TBUS]- 



Let {(wi,^i), (w2,^2)} denote the eigen decomposition of J and assume that 
^2 ^ 'Ci- It is easily seen that and .^2 respectively correspond to the maximum 
and minimum values reached by Qv when the unit vector v varies. Eigenvec- 
tors indicate the directions for which these extrema are reached. Thus, they 
respectively correspond to the average direction of image gradients over the in- 
vestigated neighborhood and to the orthogonal direction. The eigenvalues ^1 
and ^2 store information on the local coherence of the gradient field in the 
neighborhood. When (■) is implemented as a convolution with a Gaussian 
kernel Qp {p corresponds the size of the investigated neighborhood), the eigen- 
values can be combined as follows to build the required estimators s and k: 



s2~ei+6 = Tr(J) = (?,*(||VI| 



6 



(13) 



The estimator s is approximately equivalent to the average norm of the gra- 
dient. The curvature estimator is based on a second order Taylor expansion 
of / along a contour. With this approximation the eigenvalues of the struc- 
ture tensor can be expressed as functions of the strength and the curvature of 
contours. The curvature k, is then easily extracted (see Appendices for more 
details). 
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4 Experiments 



4-1 Quality of the curvature estimator 

This section illustrates the accuracy and robustness of our proposed curvature 
estimator. We investigate the influence of the sizes of the Gaussian kernels used 
to compute the gradient structure tensor and we compare our estimator with 
previous works. 

For this purpose, we generate images of ellipses with known curvature. These 
images are corrupted with different levels of Gaussian noise (see Fig. Hj). Then 
curvature is computed along ellipses with our estimator and results are com- 
pared with the true curvature. For a given noise level, the experiment is re- 
peated 40 times. The presented curves show the averages and the standard de- 
viations of estimated curvatures over this set of 40 test images. Noise levels are 
expressed using the peak signal to noise ratio defined as PSNR = 10 log 
where /max is the maximum amplitude of the input signal and a is the standard 
deviation of the noise. 

Fig. [5] illustrates the influence of the parameter a. As expected, it must be 
chosen in accordance with the noise level in the image. If a is too small, 
the direction of the image gradient changes rapidly in a neighborhood of the 
considered point. As a consequence, the second eigenvalue of the structure 
tensor increases. This explains why curvature is overestimated. 

The dependency of our estimator on the radius p of the local integration is 
depicted on Fig. [61 The presented curves show that this parameter has an 
influence only for images with strong noise. Indeed, contour information has 



18 



to be integrated over much larger neighborhoods to mitigate the influence of 
noise. 

In addition, our estimator is compared with two methods which both involve 
the computation of the second derivatives of the input image: 

• the naive operator which simply computes the curvature of isocontours of 
the image as 



• the more elaborate estimator proposed by Rieger et al. |20j . 

The latter method consists in computing the derivative of the contour ori- 
entation in the direction of the contour. Contour orientation is computed 
(modulo vr) as the eigenvector Wi of the gradient structure tensor (which 
corresponds to the largest eigenvalue). Since orientation is only known mod- 
ulo TT, the vector wi is converted into an appropriate continuous representation 
(using Knutsson mapping) prior to differentiation. 

As shown on Fig. [7] all these estimators provide fairly equivalent results along 
a contour. Note however that the naive estimator is much more sensitive to 
noise than the others. 

These estimators were also tested in places without image features. As depicted 
on Fig. [HI, both the naive estimator and the one of Rieger et al. become un- 
stable. The naive estimator fails because the denominator in falls to zero 
and because second derivatives are very sensitive to noise. Rieger's method 
can neither be used. Indeed, in a region without significant contour, the eigen- 
vector Wi of the gradient structure tensor is only determined by noise and 




(14) 
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thus exhibits rapid variations. Computing its derivatives results in a spurious 
evaluation of the curvature. This justifies the use of our estimator, which, in 
addition, requires less computations than Rieger's method since it estimates 
curvature directly from the eigendecomposition of the gradient structure ten- 
sor and does not involve their derivatives. 



4-2 Parameter selection for a new image segmentation 

Given a new image, we have to adjust some parameters to exploit at best the 
potentiahes of the proposed approach. We follow the steps below: 

(1) The image structure tensor is computed: it provides the contour inten- 
sities s, the curvatures k and the local metrics; the maximal curvature 
kmax as well as the maximal intensity Smax follow immediately. 

(2) The user chooses the minimal and the maximal edge lengths Zmin and imax 
for the model. Typically, the length Zmin is half the size of a pixel (a better 
precision has no sense given the input data) and the length /max is about 
50 times Zmin for real data. 

(3) The user then selects the reference contour intensity Sref which corre- 
sponds to reliable contours. A simple way is to visualize the thresholding 
of the image s by the value Sref, and to tune this parameter accordingly. 
It can also be automated for certain images as a percentage of Smax (typ- 
ically 90%). 

(4) After that, the procedure is the same as for classical snakes: initialisation, 
selection of energy/force parameters, evolution until rest position. 
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From the preceding paragraphs, it is clear that the proposed approach does 
not induce significantly more interaction compared with classical snakes. 

4.3 Behaviour of the deformahle model 

Adaptive vertex density is illustrated in Fig. [91 In this experiment images 
of circles with known radii are generated (left part of the figure). For each 
circle, the image space is equipped with a metric which is built as explained in 
Sect. 131 In this example Kref = 1 and Sref is computed from the input image as 
the maximum value of s over the image. Our deformable model is then used 
to segment and reconstruct the circles. Once it has converged, the Euclidean 
lengths of its edges are computed. The results are presented on the curve in 
the right part of the figure. They correspond to the expected behavior (see 

Fig. ED. 

Adaptive vertex density is also visible in Fig. [TTIfT^ As expected, changing the 
metric increases vertex density along highly curved parts of image components. 
As a result, the description of the shape of objects is enhanced while the 
number of vertices is optimized. 

Independence with respect to the resolution of input images is shown on 
Fig. [TOl Our model was tested on images of objects sampled at different rates 
(see Fig. [TTl) . As expected, the number of vertices is kept independent from 
the resolution of the input image, as far as the sampling rate ensures a proper 
representation of the highest frequencies present in the signal. If this condition 
is not satisfied, as on Fig. [121 the model uses only the available information. 
If the resolution is increased, the length of the edges of the model remain 
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unchanged, except in parts where the finer samphng rate of the image allows 
to recover finer features. 

Fig. [11] and Fig. [13] demonstrate the ability of our model to dynamically and 
automatically adapt its topology. Note that the proposed way to build the met- 
ric is especially well suited to objects with thin and elongated features. With 
previous approaches p][2] automated topology changes can only be achieved 
using grids with a uniform resolution determined by the thinest part of objects. 
Their complexities are thus determined by the size of the smallest features to 
be reconstructed. The involved computational effort is therefore wasteful since 
much more vertices are used than required for the accurate description of ob- 
jects. In contrast, replacing the Euclidean metric with a metric designed as 
described in this paper virtually broaden thin structures. As a consequence, 
even for large values of 6, the inequality dSj) where cIe has been replaced by dji 
is not satisfied for two vertices u and v located on opposite sides of long-limbed 
parts of image components. Self-collisions are thus detected only where they 
really occur. At the same time, the number of vertices is kept independent 
from the size of the finest details to be recovered. 

Fig. [T31[T^ illustrate the behavior of our deformable model on biomedical im- 
ages. The input image (Fig. [13] top-right) is a fiuorescein angiogram that re- 
veals the structure of the circulatory system of the back of an eye. In addition 
to the classical regularizing forces, the vertices of the active contour undergo an 
application-specific force designed to help recovering blood vessels. This force 
pushes vertices in the direction of the outer normal and stops when the local 
gray level of the image falls under the average gray level over a neighborhood. 
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More formally, the force undergone by a vertex v is defined as 



F, = {I{v)-{gr*I){v))xn^ 



(15) 



where / is the input image, gr is a Gaussian filter used to estimate the average 
gray level over a neighborhood, and is the outer normal to the deformable 
curve at vertex v. The presented results demonstrate the possibility to use 
additional forces designed to extend or improve deformable models [T3lll2fl4j . 




Furthermore, computation times are given in Table [2J They show that re- 
ducing the number of required iterations and the number of vertices largely 
compensate for the time used to compute of the metric. These computation 
times are given in Table [3] for different size of input images. 

At last, since the space is expanded only in the vicinity of image contours, 
vertices travel faster in parts of the image without feature. When approaching 
object boundaries, the deformable curve propagates slower and progressively 
refines. The advantage is twofold. First, the number of iterations required for 
the model to reach its rest position is reduced. Second, the cost of an iteration 
is reduced for parts of the deformable curve far away from image features, 
namely parts with a low vertex density. By this way a lot of computational 
complexity is saved when the deformable model is poorly initialized. This is 
especially visible on Fig. [11] (right) where the position of the model has been 
drawn every 50 iterations. 

5 Conclusion 

We presented a deformable model that adapts its resolution according to the 
geometrical complexity of image features. It is therefore able to recover finest 
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Table 2 

Number of iterations and time required to reach convergence. The table also indi- 
cates how many vertices are used to represent the whole shapes shown in Fig. [TH 
The two last columns describe the minimum and maximum length (in pixels) of an 



ige of the c 


eformable n 
iterations 


lodel. 
total time (s) 


vertices 


min. edge 
length 


max. edge 
length 


uniform 
adaptive 


350 
250 


78.5 
37.35 (+1.24) 


3656 
2065 


0.5 
0.35 


1.25 

25 



Table 3 

Computation times required to build the metric for different sizes of input images. 



resolution of input image 


100 150 200 


250 300 350 


400 


computation time (s) 


0.16 0.36 0.65 


1.00 1.45 1.98 


2.58 



details in images with a complexity almost independent from the size of input 
data. Admittedly, a preprocessing step is required to build the metric. How- 
ever, involved computational costs are negligible and, as a byproduct, these 
precomputations provide a robust gradient estimator which can be used as a 
potential field for the deformable model. Most of the material used in our pre- 
sented deformable model has a straightforward extension to higher dimensions 
[T6f2T] . We are currently working on the extension of the whole approach to 
3D images. 

A Second order approximation of contours 

In this section, we consider a contour that is tangent to the at the 

origin. This is expressed as lx{0,0) = and ly{0,0) = s, where Ix, ly and s 
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denote the partial derivatives of / and the strength of the contour. 



Prom the definition of a contour as a maximum of the norm of the gradient 

(0,0) 



in the gradient direction follows ^ = 0. Once expanded, this leads to 



7,,(0,0) = 0. 

Let t and n denote the vectors tangent and normal to the investigated contour: 
n = n^Yii and t = n^. From the definition of curvature follows ^ = Kt. 
Replacing t and n by their expression as functions of /, and then expanding 
and evaluating this expression at point (0, 0) yields Ixxi^, 0) = sk. 

From the above statements we get a second order Taylor expansion of / as 
I{x, y) - /(O, Q)^s(y + ^^Kx"^ + hyxy + o(x^ y^) . (A.l) 



In addition, if we assume that the strength of the contour remains constant 

(0,0) 



along the contour, we get ^^q^^P =0. Expanding the previous expression 



leads to lxy{0, 0) — 0. 

With this additional hypothesis, / may be rewritten as 



I{x, y) - /(O, 0) = s ( 2/ + ]^KX^ ) + o{x\ y^) . (A.2) 



B Structure tensor of a peirabolic contour 

In this section we compute the eigenvalues of the structure tensor along a 
contour with strength s and with local curvature k. 
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B.l Contour with constant intensity 



Following approximation ( ]A.2I) we consider the image / defined as 

I{x,y) = s(y + ^Kx^^ . (B.l) 

For symmetry reasons, we know that the eigenvectors of the structure tensor 
J at point (0, 0) are aligned with the x axis and y axis. In this special case, 



the eigenvalues of J are given as = ly = s"^ and ^2 = Ix = s'^k^x^. If 
the averaging operation over a neighborhood is implemented as a convolution 
with a Gaussian function gp, this yields ,^1 = and ,^2 = s^/t^p^. In practice 
only ,^1, ^2 and p are known. Curvature (up to sign) is easily computed from 
these quantities as 

\k\ = -M (B.2) 
B.2 Contour with varying intensity 



In this subsection we show that the estimator described in the previous para- 
graph remains valid to estimate the curvature of a contour with a varying 
intensity 



We start with equation flA.ip : 

1 



I{x, y) = s[y + ^Kx^ ) + hyxy , (B.3) 



from which we get ^1 = + /J and ^2 = s^K^p^ + I^yP 



The curvature estimation k is thus written : 



1/^2 / K^S^ + Ixy 



26 



If K = we get 



ut) + 



(B.5) 



Assuming that the contour intensity is significantly greater than its linear 
variation along x, we obtain ^ ~ 0. Replacing in fIB.SP we get k ~ 0. 



If K 7^ 0, we get 



(B.6) 



As shown before ^ ~ for a reliable contour. As a consequence, such a 
contour k ~ k, which shows that the curvature estimator remains available 
for contours with a varying intensity. 



C Implementation issues 

The gradient structure tensor is implemented as successive convolutions of the 
input image with two Gaussian functions g^, Qp and their partial derivatives: 

J = gp*{V{I*g„)x'V{I*g„)) . (C.l) 

Convolutions are implemented efficiently as a product in the frequency domain 
and could be further improved using recursive implementations of Gaussian 
filters 

The parameter a determines how much the image gets smoothed before com- 
puting its derivatives. It is thus chosen in accordance with the noise level in the 
image. The parameter p determines the size of the neighborhood over which 
the gradient information is integrated. The infiuences of these parameters are 
studied experimentally in Sect. 14. 1[ 
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Since the metric has to be computed everywhere in the image our estimator 
must remain stable in regions without contours {i.e. in regions where ~ 0. 
Therefore k is computed as follows: 



where e is an arbitrary positive constant. By this way the denominator never 
vanishes and k, falls to in places without image structure. In the vicinity of 
image contours e may be neglected in front of and we get back estimation 
( IB. 21) . Experimentally, a suitable choice for this constant is e = jq^T'^^ where 
^max f|g]2otes the maximum value of over the image. 
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Fig. 1. Illustration of the proposed approach to shape extraction. Top row: ex- 
traction with uniform sampling of the model. The edge length has approximately 
the pixel width. Computation statistics: 910 iterations, 10.14s, 458 vertices. Middle 
row: coarse to fine extraction. At each level, the edge length has approximately 
the pixel width. Computation statistics: 310+810+760+1020+510 iterations, 9.23s 
= 0.08+0.40+0.87+3.10+4.78, 392 vertices. Bottom row: extraction with adaptive 
sampling of the model. The edge length varies between half the pixel width and 
20 times the pixel width. Computation statistics: 280 iterations, 1.73=0. 31+1. 42s 
(precomputation + evolution), 150 vertices. 
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Fig. 2. Edge length (up to a factor at most C,) depending on the strength s and cur- 
vature K of the underlying contour. It is assumed that edges run along the contour. 
In region A contours are too weak or too straight. Therefore, edges keep their max- 
imum length. In region B contours are considered as reliable and have a curvature 
higher than the reference curvature Kref- The length of the edges increases linearly 
with the radius of curvature of underlying contours. In region C contours have a 
varying reliability and have a curvature higher than the reference curvature. The 
length of edges depends on both s and k. The separation between regions A and C 
corresponds to contours for which — — — = 1. It corresponds to contours for which 
curvature and/or strength fall too low to let the model increase its vertex density 
safely. 
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Fig. 3. Left: edge length (up to a C factor) as a function of the radius of curvature of 
the underlying contour. The solid line corresponds to a reliable contour (s > Srcf), 
the dashed line corresponds to a weaker contour (s < Sref)- The hatched part of 
the graph cannot be reached sincG estinicited. curvaturGS cannot cxcggcI ^max* 

Right: 

edge length (up to a factor at most Q as a function of the strength of the underlying 
contour. The solid line corresponds to a contour with the highest possible curvature 
(k = Kmax)- The dashed lines corresponds to less curved contours (k < Kmax)- In 
both figures, it is assumed that edges run along the contour. 




Fig. 4. Images used to test curvature estimators (from left to right PSNR = 40 dB, 
30 dB and 20 dB). The curvature is estimated along the border of ellipses and are 
compared with the true curvature for different estimators and different values of the 
parameters. Results are presented on Fig. [5][71 
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Fig. 5. Estimated curvature as a function of the true curvature for different values 
of a and for p = 10 (the x-axis and y-axis respectively correspond to the true 
and estimated curvatures). The three graphics correspond to different noise levels: 
PSNR = 40 dB (top-left), PSNR = 30 dB (top-right) and PSNR = 20 dB 
(bottom) . 
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Fig. 6. Estimated curvature as a function of the true curvature for a = 5 and 
for different values of p (abscissa and ordinate respectively correspond to the true 
and estimated curvatures). The three graphics correspond to different noise levels: 
PSNR = 40 dB (top-left), PSNR = 30 dB (top-right) and PSNR = 20 dB 
(bottom) . 




Fig. 7. Comparison of curvature estimators along the contour of a noisy ellipse (left 
PSNR = 40 dB, right PSNR = 20 dB). (A) our estimator, (B) Rieger's estimator, 
(C) naive estimator. 
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Fig. 8. Comparison of curvature estimators. Left: test image {PSNR = 40 dB). 
Right: curvature estimated along the radius drawn on the left figure. (A) our esti- 
mator, (B) Rieger's estimator. The solid line represents the inverse of the distance 
to the center of the circle and the arrow indicates the position of the contour. The 
naive operator is not displayed since it is too unstable. 
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Fig. 10. Final number of vertices on the curve depending on the resolution of 
input image. The segmentation/reconstruction results as well as the evolution of 
the model are shown on Fig. [TTJ 




Fig. 11. Left, center: reconstruction of identical objects sampled at resolutions 
40x40 and 100x100. Right: evolution of the deformable model every 50 iterations. 
The outer square corresponds to the initial position of the model. 
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Fig. 12. Segmentation/Reconstruction of the same object sampled at increasing 
resolutions (50x50, 100x100, 200x200 and 400x400). For the four images all the 
parameters used to build the metric or attract the model toward object boundaries 
are identical. Please note that the deformable model automatically adapts its reso- 
lution to represent available image features as well as possible, while optimizing the 
number of vertices. 
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Fig. 13. Segmentation process on an angiography. Top-left: input image. Other 
images: steps of the evolution of the deformable curve. The model is driven by 
an inflation force which stops when the local gray level is lower that the average 
gray-level in a neighborhood. Please note the topology changes when parts of the 
deformable curve collide. Computation times are given on Table [2] 



39 




Fig. 14. Segmentation of the angiography shown on Fig. [T3j Left: results without 
adaptation. Right: results with a metric built as described in [3l Top: final result. 
Bottom: detailed view in a region which exhibits much adaptation of the vertex 
density as well as a complex topology. Please note how the length of edges is adapted 
according to the structures found in the input image and how this enables the 
deformable model to enter small gaps and recover structures finer than its edges. 
Computation times are given in Table [2j 
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